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(57) Method of automatic registration of two- and 
three-dimensional angiography images by comparison 
of a two-dimensional digital subtracted Oangiography 
image with data on a three-dimensional image recon- 
structed from rotational angiography sequences, in 
which a field of distortions in the image is estimated, a 
conical projection matrix is estimated and an approxi- 
mation is made of a rigid transformation in space equal 
to the difference between an initial registration based on 
the field of distortions and on the conical projection 
matrix and a perfect registration. 
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Description 

[0001] The present invention concerns the field of image processing, particularly of two- and three-dimensional 
radiological images. 

5 [0002] As is known, radiological apparatuses comprise a means of emission such as an X-ray tube and a means of 
reception of the emission, such as a solid state detector or even a scintillator and a video camera, of CCD type, for 
example. 

[0003] The means of emission and the means of reception of X-rays are generally supported by a mobile system 
with one or more axes, making it possible to take pictures at different angles of incidence. The means of reception is 
w connected to image processing means making possible the generation of three-dimensional images from a series of 
two-dimensional images picked up by the means of reception. 

[0004] It is important to obtain a good match between a three-dimensional image and a two-dimensional image 
taken upon a stage in which the patient undergoes a particular procedure such as placement of a catheter in the field 
of angiography, in order to be able to follow the movement of the catheter in the two-dimensional image, but also in the 

75 three-dimensional image. 

[0005] A method of automatic registration of subtracted angiography images is known from the article "Fully auto- 
matic 3D/2D Subtracted Angiography Registration- by KERRIEN, LAUNAY, BERGER, MAURINCOMME, VAILLANT 
and PICARD, published in CARS'99, Proceedings of the 13 lh International Congress and Exhibition, Paris, June 23-26, 
1999. This article indicates that reconstructed three-dimensional images of cerebral vascularization can be obtained 

20 from a large number of subtracted angiography images and that in interventional neurology it is worthwhile to register 
the three-dimensional volume reconstructed on the instantaneous two-dimensional images. 

[0006] In interventional radiology, it is desirable for the practitioner to be able to know at any time that the catheter 
is in the patient's body with the utmost possible precision. The practitioner can at the present time deduce this informa- 
tion from digital subtracted angiography images called DSA, which he associates with preoperative magnetic reso- 

25 nance thanks to his anatomical knowledge. DSA images offer, among other things, real time imaging with very high 
space resolution. Three-dimensional reconstructions of blood vessels called "3DXA" have also been used recently from 
rotational angiography sequences made by rapid rotation of the X-ray tube and of the camera over half a turn and the 
taking of about fifty DSA images, which are the projections on input of a tomography algorithm producing the 3DXA 
image on output. For more information on this technique, the reader is invited to refer to the thesis of LAUNAY, "Local- 

30 ization and 3D reconstruction from stereotaxic angiograms," doctoral thesis, National Polytechnic Institute of Lorraine, 
Nancy, France, 1996. 

[0007] These reconstructions make possible a very good appreciation of angioarchitecture. Furthermore, those 
three-dimensional images can be used in real time according to several types of visualization, such as maximum inten- 
sity projection, isosurface, volume melting, virtual endoscopy or even reformatted cross -section, and are a further assist 
35 to the diagnoses of practitioners. 

[0008] The present invention concerns an improved method of registration. 

[0009] The present invention also concerns a method of registration of millimetre or submillimetric precision at 
short calculation time. 

[0010] The method of automatic registration is intended for two- and three-dimensional angiography images by 
40 comparison of a two-dimensional digital subtracted angiography image with data on a three-dimensional image recon- 
structed from rotational angiography sequences, in which a field of distortions in the image is estimated, a conical pro- 
jection matrix is estimated and an approximation is made of a rigid transformation in space equal to the difference 
between an initial registration based on the field of distortions and on the conical projection matrix and a perfect regis- 
tration. With the rigid transformation comprising a translation part and a rotation part, an approximation is made of the 
45 translation part by considering the rotation part as known, the optimal transformation being attained when the correla- 
tion of the digital subtracted angiography image and the reconstructed three-dimensional image data is at a maximum. 
[001 1 ] It is assumed preferably that the residual registration error is due to an error of positioning on translation and 
rotation of the reconstructed three-dimensional image in the mark of an imaging means. 

[001 2] Preferably, the difference between the reconstructed three-dimensional image data and the two-dimensional 
so digital subtracted angiography image is considered due to a slight displacement following a rigid transformation, the 
rigid transformation being determined by a modified optical flux technique. 

[0013] The reconstructed three-dimensional image data preferably relate to a two-dimensional image obtained by 
conical projection of maximum intensity of the reconstructed three-dimensional image. This type of projection links a 
single voxel to each pixel. The depth of the voxel chosen as value of one of the space coordinates at a given instant can 
55 be used for each pixel on projection. 

[0014] Advantageously, on each displacement of the digital subtracted angiography image, the new two-dimen- 
sional image is calculated from the old two-dimensional image, without new calculation of projection of the recon- 
structed three-dimensional image. The duration of calculation necessary is considerably reduced from several hours to 
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one or two minutes. 

[0015] In an embodiment of the invention, the initial registration is given by a marker unit of known three-dimen- 
sional position. 

[0016] In another embodiment of the invention, the initial registration is given by calibration of the angiography 
5 machine. 

[0017] Thus, starting from an initial registration, the initial error in translation is reduced by maximization of the cor- 
relation between the digital subtracted angiography image and the reconstructed three-dimensional image data, and by 
refinement of registration by a modified optical flux technique, it is possible to obtain a registration of millimetric or 
inframillimetric precision perfectly sufficient for real use in the course of an operation, notably, of interventional neuro- 
w surgery or radiology, at the end of a calculation time limited to a few minutes, and even to a few tens of seconds. 

[0018] The invention applies, in particular, to radiology with three-dimensional images in which the blood vessels 
are visible. The informational content of such images is slight. 

[001 9] The present invention will be better understood by studying the detailed description of an embodiment given 
by way of n on limitative example. 
is [0020] The invention makes possible a matching of DSA images and 3DXA images. 

[0021] "DSA image" means here the image of maximum opacification up to row N in the acquired sequence, that 
is, each pixel of the resultant image takes the smallest value encountered on the N first images of the sequence, or the 
image of row N in the acquired sequence. Row N of the Image is either chosen by the user or fixed in relation to the rate 
of acquisition. 

20 [0022] Furthermore, the 3DXA volume is projected in view called conical "MIP" (maximum intensity projection). 
Each voxel of the 3DXA volume can be considered a sphere of radius equal to the size of the voxel, and is projected in 
the image plane. The final value of each pixel is the maximum value projected in that place. The image thus obtained 
is described as the "MIP image." As a variant, it is possible to use a conical projection by addition of the values of voxels 
in each pixel. 

25 [0023] Registration between the two image is made by estimate of a field of distortions of the image, a technique 
known from the article "Machine precision assessment for 3 D/2D digital subtracted angiography images registration' by 
KERRIEN, VAILLANT, LAUNAY, BERGER, MAURINCOMME and PICARD, in SPIE Medical Imaging, Volume 3338, 
pages 39-49, 1998. 

[0024] After estimate of the field of distortions in the image, an estimate is made of a conical projection matrix. A 
30 calibration makes it possible to find out the intrinsic parameters of this matrix, the extrinsic parameters being deter- 
mined as explained below. 

[0025] The radiology machine, once calibrated, supplies an initial registration which differs from the registration 
sought of perfect registration by a rigid transformation (rotation + translation) in three-dimensional space. The rigid 
transformation possesses a rotation part of low amplitude, while the translation part can be considered unknown. 

35 [0026] Initially, the translation part is found by considering the rotation part of the projection matrix perfectly known. 
The optimal position is reached when the centered standardized correlation of the DSA image with the MIP image is 
maximal. It is assumed that the residual registration error is due to a slight positioning error in both rotation and trans- 
lation of the 3DXA volume in the camera mark corresponding to the registration obtained. The MIP image is then con- 
sidered taken at a time t 1t since the 3DXA volume moves according to a slight rigid transformation. An image is 

40 acquired at time t 2 and the DSA image is obtained. The apparent movement observed between the present MIP image 
and the DSA image gives information concerning the slight rigid transformation sought. This information is extracted by 
a modified optical flux technique. 

[0027] The centered standardize correlation between the DSA image ! d , and the MIP image, l m , is defined by: 



so where D is the domain common to both images and ll d , ! m are the means of images l d and l m respectfully. The current 
registration is then considered given by the projection matrix M. A new projection matrix M' is sought, which differs from 
M by the translation D = (D x , Dy, D z ) . Considering a point P which is projected in the coordinate locus (u, v) in the 
initial image obtained by projection of the volume according to the matrix M and in the coordinate locus (u\ v') in the 
final image obtained after translation, we then have: 



45 
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M'P = 




(I) 



[0028] The matrix M is then standardized by Toscani's method ("System of calibration and perception of movement 
in artificial vision," G. Toscani, doctoral thesis, Universite de Paris Sud, Orsay, 19B7) and one deduces therefrom: 



[0029] This equation makes possible an interpretation of translation D. D x and D y produce a translation of the initial 
image, while the term D 2 creates a zoom. This homothety Is centered on the upper left-hand corner of the image and 
is accompanied then by a translation of structures in the image. The dependence on parameters can be detrimental to 



[0030] If, now, we have a translation (du, dv) parallel to the image plane associated with a magnification relative to 
the center of the image (u c v c ), we obtain: 



[0031] Those parameters are independent and lend themselves well to optimization of the correlation thanks to the 
fact that they are more intuitive. They are not, however, equivalent to the vector D. In fact, s depends on the coordinates 
X, Y and Z of point P. We therefore considered that s was constant throughout the volume, which means considering 
that this volume is punctual. It can be estimated that reconstruction of the 3DXA volume occupies a sphere, for exam- 
ple, 15 cm in diameter, situated midway between the focal point and the image plane, the focal distance of which is, for 
example, in the order of one meter. The approximation is sufficient, for the matrices obtained indicate that the variation 
of s is in the order of 1 %. The hypothesis is therefore valid and the two equations above are identified in order to deter- 
mine a bijection between the vector D and the triplet (du, dv G). 
[0032] The procedure of optimization of these three parameters is as follows: 

- exhaustive low-resolution search, for example, 64x64 pixels. The limits of variation of the parameters are either set 
by the angiography machine or respond to a broad criterion, which is that the two images overlap over at least one- 
quarter of the surface; 

- maximum resolution search, which consists of an exhaustive search along the axes of the mark of space of the 
parameters. Only one parameter at a time is modified to generate the search space. In the case of three parame- 
ters, we therefore have only six correlation calculations instead of 26 to be made. It is to be noted, finally, that, 
parameter G being very different from du and dv in amplitude and influence, we alternate between its optimization 
and the common optimization of du and dv. The parameters acting in the axial plane (with high resolution) and 
those exerting their influence in the orthogonal direction (with lower resolution) are thus distinguished. 

[0033] The parameters du, dv and G, which maximize the correlation between the DSA image and the MIP image, 
being found, the projection matrix M is modified accordingly and the 3DXA volume is projected to generate a new MIP 
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image. That new image is close enough to the DSA to be able to lend itself to an optical flux calculation. The difference 
between the MIP image and the DSA image is explained by an error on the position of the 3DXA volume in the mark of 
the camera. This error is akin to a slight rigid displacement. 

[0034] Let us consider a point P = (X, Y, Z) in space at time t and point P' = (X\ Y\ 2') reached at time t* from P 
5 by a slight rigid displacement, consisting of rotation R = R A R B R C (A, B and C being the angles of rotation around the 
base vectors of the three-dimensional mark) and translation T = (U, V, W) . We then have: 



P' = RP + T 



io [0035] In the hypothesis where the movement is slight, we have: 



P'-P = P = ftXP + T with Q = (A, B, C) 



[0036] In the internal mark of the camera, point P is projected on the pixel of coordinate (u, v) of the image accord- 
is ing to the formulas: 

X Y 
u ~ a ^ and v = a ^ 

where a is the ratio of the focal distance to pixel size for the camera. By deriving that equation in relation to time and 
20 combining it with the previous one, we obtain: 



25 



u = qlB - Cv uv + — + 

a a Z 

v = - qA + Cu+ *uv - ±v> + aV ' Wv 
a a Z 



(6) 



30 

[0037] This equation corresponds almost to a transformation in the image plane, except for the term in 1/Z. This is 
the problem raised by resolution of the three-dimensional movement from the apparent movement of the equation defin- 
ing P" - R Only five parameters out of the six can be resolved, for an uncertainty remains as to the depth of the object 
in the field given by the variable Z. An object in translation parallel to the image plane will generate the same apparent 

35 movement as an object twice as small, situated twice as close to the optical center and driven by a movement parallel 
to the first, but twice as weak. We can, however, take advantage of the projection of MIP type. In fact, that projection 
joins a single voxel to each pixel. For each pixel, on calculation of the MIP projection, we can, consequently, join the 
depth of the voxel chosen as value of Z. The equation determining the derivatives of u and v then completely links the 
coordinates of the pixels in the image to the six parameters of the movement sought. 

40 [0038] The values of u and v for each pixel remain to be determined. According to the hypothesis of optical flux, the 
intensity of the object between the DSA image and the MIP image is constant. That is not absolutely observed, in the 
sense that the two images are obtained by different means and do not represent exactly the same objects, a real object 
in one case and a tomographic reconstruction in the other. Following this idea, however, one writes: 

S-»'-CH- 



[0039] By combining that equation with the previous equation, we obtain for each pixel two equations that must 
so respect the six parameters of the movement. We then arrive at a system widely overdetermined with six unknowns, 
which are resolved to the least squares by determination of the pseudo-inverse. 

[0040] The two phases explained above are based on a large number of MIP projections of the volume. This is evi- 
dent in the correlation optimization procedure. 

[0041] In the case of the modified optical flux, the results are qualitatively good, but quantitatively bad, the ampli- 
55 tude of the movement being very much underestimated, even though its direction may be correct. An iterative resolution 
of the slight residual movement is then used, which requires a number of MIP type projections, namely, one per itera- 
tion. 

[0042] Now, the generation of a conical MIP type projection of a 3DXA volume takes in the order of one second. For 
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example, in case the low-resolution image is ot size 64*64, the number of projections required on exhaustive calculation 
ot the low-resolution correlation is the product of the number of values of parameters du, dv and G, respectively equal 
to 32, 32 and 35. Therefore, 35840 correlation calculations must be made and, consequently, just as many MIP type 
projections. At the rate of one second per projection calculation, the duration becomes excessive. 
5 [0043] However, the equation defining u' and V and that defining u and v each define a transformation in the image 
plane, if we know the parameters of the associated movement (du, dv and G for the correlation; A, B, C, U, V and W for 
the optical flux). Thus, considering a movement of translation or slight rigid displacement, we can calculate the new MIP 
image from the old one, without going through a reprojection of the 3DXA volume. The calculation time is considerably 
improved. 

w [0044] In the case of modified optical flux, the transformation in the image plane cannot be reversed because of the 
term in 1/Z. It follows that the images undergo an artifact. The degraded zones, that is, the pixels not reached by the 
transformation, can, however, be localized and they can be avoided in establishment of the system of equations. 
[0045] The invention therefore makes possible an entirely automatic registration between the DSA image and the 
3DSXA volume. A test on an UltraSparc® station rated at 200 MHz was carried out in 1 .5 minutes. The reference reg- 

15 istration can be given by using a stereotaxic frame. Considering two DSA images taken roughly perpendicular and reg- 
istered with a 3DXA volume, particular structures can be targeted, such as branches of arteries, pronounced curves, 
etc., in order to reconstruct the position of the cursor so defined in the 3DXA volume. 

[0046] In other words, the registration of DSA images with volume images of blood vessels is carried out by calibra- 
tion of the angiography machine or by use of markers, reduction of the initial error by estimate of a three-dimensional 
20 translation, by maximization of the correlation between the DSA image and the conical projection image of the volume, 
MIP, for example, and use of a modified optical flux technique to refine the registration. 

[0047] It is therefore possible to precisely locate a catheter in a 3DXA volume, to visualize a blood flow in the 3DXA 
volume, to estimate the patient's movement in relation to the angiography machine during an image acquisition and to 
eliminate use of the stereotaxic frame on treatment of a patient. It is also possible to avoid the use of presegmentation 
25 of images, which is difficult to carry out and consumes large calculation capacities and requires a manual stage of ver- 
ification of results by the user. Furthermore, segmentation disrupts the results of registration because of errors, super- 
positions of arteries in the DSA image, bifurcations, tangencies of arteries, etc., and which can then be detrimental to 
precision. A high degree of precision is attained by using variables calculated from the base data of the image. 

30 Claims 

1. Method of automatic registration of two- and three-dimensional angiography images by comparison of a two- 
dimensional digital subtracted angiography image with data on a three-dimensional image reconstructed from rota- 
tional angiography sequences, in which afield of distortions in the image is estimated, a conical projection matrix 

35 is estimated and an approximation is made of a rigid transformation in space equal to the difference between an 
initial registration based on the field of distortions and on the conical projection matrix and a perfect registration, 
the rigid transformation comprising a translation part and a rotation part, an approximation is made of the transla- 
tion part by considering the rotation part as known, the optimal transformation being attained when the correlation 
of the digital subtracted angiography image and the reconstructed three-dimensional image data is at a maximum. 

40 

2. Method according to claim 1 , in which it is assumed that the residual registration error is due to a positioning error 
on translation and on rotation of the reconstructed three-dimensional image in the mark of an imaging means. 

3. Method according to claim 2, in which the difference between the said data on the reconstructed three-dimensional 
45 image and the two-dimensional digital subtracted angiography image is considered due to a slight displacement fol- 
lowing a rigid transformation, and said rigid transformation is determined by a modified optical flux technique. 

4. Method according to any one of the foregoing claims, in which the said data on the reconstructed three-dimensional 
image consist of a two-dimensional image obtained by conical projection of maximum intensity or by addition of the 

50 values of voxels of the three-dimensional reconstructed image. 

5. Method according to claim 4, in which, on each displacement of the digital subtracted angiography image, the new 
two-dimensional image is calculated from the old two-dimensional image, without any new calculation of projection 
of the reconstructed three-dimensional image. 

55 

6. Method according to any one of the foregoing claims, in which the initial registration is given by a marker unit of 
known three-dimensional position. 
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7. Method according to any one of claims 1 to 5, in which the initial registration Is given by calibration of the angiog- 
raphy machine. 
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